/*
 * kmeans.c
 *
 *  Created on: 3 nov. 2012
 *      Author: Rania ToumiGasri
 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include "kmeans.h"
#include "utils.h"
#include "configuration.h"

Map* createStructureOfCentres(int dimensionDesDonnees) {
	Map* mapCentre;
	mapCentre = (Map*) malloc(sizeof(Map));
	mapCentre->nombreDeCentre = config.centerNumber;
	mapCentre->dimensionDesDonnees = dimensionDesDonnees;
	// allocation dynamique du tableau ayant les centres.
	double** centreTab = (double**) malloc(
			config.centerNumber * sizeof(double*));
	int i, j;
	for (i = 0; i < config.centerNumber; i++) {
		//la derniere case du tableau contient le nombre de points de la classe i.
		centreTab[i] = (double*) malloc(
				(dimensionDesDonnees + 1) * sizeof(double));
		for (j = 0; j <= dimensionDesDonnees; j++) {
			centreTab[i][j] = 0;
		}
	}
	mapCentre->codeBook = centreTab;
	return mapCentre;
}
////////////////////////////////////////////////////:
Map* initialiserLesCentresMobiles(Data* myData) {
	Map* mapCentre = createStructureOfCentres(myData->dimensionDesDonnees);
	// initialiser le tableau des matrices de covariances
	covariances = (double***) malloc(config.centerNumber * sizeof(double *));
	int k, l, m;
	for (k = 0; k < config.centerNumber; k++) {
		covariances[k] = (double**) malloc(
				myData->dimensionDesDonnees * sizeof(double *));
		for (l = 0; l < myData->dimensionDesDonnees; l++) {
			covariances[k][l] = (double*) malloc(
					myData->dimensionDesDonnees * sizeof(double));
			for (m = 0; m < myData->dimensionDesDonnees; m++) {
				covariances[k][l][m] = 0;
			}
		}
	}
	int i, nb, j;
	// prendre aléatoirement des points de l'ensemble de points dans Data
	/*tab contient les indices des centres qu'on  va tirer aléatoirement entre 0 et nbData*/
	int tab[config.centerNumber];	// contient les indices des centres
	for (i = 0; i < config.centerNumber; i++) {
		do {
			nb = rand() % myData->nombreDesDonnees;
			tab[i] = nb;
		} while (existingCenter(tab, nb, i));
		tab[i] = nb;
	}
	// prendre les donnees proprement dite
	for (i = 0; i < config.centerNumber; i++) {
		nb = tab[i];
		for (j = 0; j < myData->dimensionDesDonnees; j++) {
			mapCentre->codeBook[i][j] = myData->data[nb][j];
		}
	}
	return mapCentre;
}
Map* apprendreLesCentresMobiles(Map* myMap, Data* myData) {
	int i = 0;
	Map* myNewMap = myMap;
	Map* myOldMap;
	do {
		myOldMap = myNewMap;
		attribuerLesDonneesAuxCentresMobiles(myOldMap, myData);
		myNewMap = calculerLesCentresDeClasses(myOldMap, myData, NULL );
		i++;
	} while (!convergenceDetecte(i, myOldMap, myNewMap));
	return myNewMap;
}

void attribuerLesDonneesAuxCentresMobiles(Map* myMap, Data* myData) {
	int i;
	for (i = 0; i < myData->nombreDesDonnees; i++) {
		myData->data[i][myData->dimensionDesDonnees] = (double) nearestCenter(
				myData->data[i], myMap);
	}
}
/**
 * calcule les coordonnées des nouveaux centres et renvoie
 * une nouvelle Map* qui contient les centres mis à jour.
 */
Map * calculerLesCentresDeClasses(Map* myMap, Data* myData, int* bmus) {
	//creer une nouvelle structure qui contientdra les nouveaux centres.
	Map* myNewMap = createStructureOfCentres(myMap->dimensionDesDonnees);

	int i, j;

	//pour chaque point
	for (i = 0; i < myData->nombreDesDonnees; i++) {
		//le point  (i.e. tableau des coordonnées du point)
		double * point = myData->data[i];
		//l'index de la classe (centre mobile) à laquelle le point appartient
		int indexCenter = point[myData->dimensionDesDonnees];
		//le centre dans la nouvelle structure
		double * centre = myNewMap->codeBook[indexCenter];
		//le nombre de points de la classe
		double nbPoints = centre[myMap->dimensionDesDonnees];
		//pour chaque coordoonée du point calculer la coordonnée du noueau centre.
		for (j = 0; j < myMap->dimensionDesDonnees; j++) {
			//coordCentre = sum(coorpoints)/nbpoint
			centre[j] = (point[j] + (centre[j] * nbPoints)) / (nbPoints + 1);
		}
		nbPoints++;
		centre[myMap->dimensionDesDonnees] = nbPoints;
	}
	return myNewMap;
}
/**
 * la convergence est detecté si le nb d'iteration depasse le nombre maximale
 *  ou si l'affectation des points ne change plus.
 */
int convergenceDetecte(int nbIteration, Map* oldCentres, Map* newCentres) {
	int result = 1;
	if (nbIteration <= config.maxIteration) {
		int i, j;
		for (i = 0; i < oldCentres->nombreDeCentre; i++) {
			for (j = 0; j < oldCentres->dimensionDesDonnees; j++) {
				if (oldCentres->codeBook[i][j] != newCentres->codeBook[i][j]) {
					result = 0;
					break;
				}
			}
		}
	}
	if (!result) {
		printCenters(oldCentres);
		printCenters(newCentres);
	}
	return result;
}
/*retourne le centre le plus proche*/
int nearestCenter(double* point, Map* myMap) {
	int i, j, x;
	double distance;
	double tab[config.centerNumber];
	for (i = 0; i < config.centerNumber; i++) {
		distance = 0;
		double * centre = myMap->codeBook[i];
		for (j = 0; j < myMap->dimensionDesDonnees; j++) {
			distance += ((point[j] - centre[j]) * (point[j] - centre[j]));
			printf(" distance = %f = (%f-%f)*(%f-%f)\n", distance, point[j],
					centre[j], point[j], centre[j]);
		}
		tab[i] = sqrt(distance);
		printf("la distance entre ");
		printPoint(point, myMap->dimensionDesDonnees);
		printf(" et le centre ");
		printPoint(centre, myMap->dimensionDesDonnees);
		printf(" egale %f\n", distance);

	}
	x = minimumDistance(tab, config.centerNumber);
	char centre[256];
	sprintf(centre, "(%f,%f)", point[0], point[1]);
	printf("le centre le plus proche de %s est le centre numero %d\n", centre,
			x);
	return x;
}
// Cette méthode retourne la classe ayant la distance minimum (l'indice du tableau est le nombre de classe).
int minimumDistance(double* T, int n) {
	int i, min = 0;
	for (i = 1; i < n; i++) {
		if (T[i] < T[min]) {
			min = i;
		}
	}
	return min;
}

//calculer l'ecart type de chaque classe.
void calculateCovariances(Map* centres, Data* donnees) {
	int i, j, k, l;
	double * centreCourant;
	double * pointCourant;
	for (i = 0; i < config.centerNumber; i++) {
		centreCourant = centres->codeBook[i];
		for (j = 0; j < donnees->nombreDesDonnees; j++) {
			pointCourant = donnees->data[j];
			if (pointCourant[donnees->dimensionDesDonnees] == i) {
				for (k = 0; k < donnees->dimensionDesDonnees; k++) {
					for (l = k; l < donnees->dimensionDesDonnees; l++) {
						double produit = (pointCourant[k] - centreCourant[k])
								* (pointCourant[l] - centreCourant[l]);
						covariances[i][k][l] += (double) (produit
								/ centreCourant[centres->dimensionDesDonnees]);
					}
				}
			}
		}
	}
}

void generateGnuPLotData(Map* centres, Data* donnees) {
	FILE *centremobilefile = fopen("centremobiles.txt", "w");
	printCentersToFile(centres, centremobilefile);
	fclose(centremobilefile);
	FILE *covariancesFile = fopen("covariances.txt", "w");
	int i, j, k;
	double ** matriceCovariance;
	for (i = 0; i < config.centerNumber; i++) {
		fprintf(covariancesFile, "classe %d :\n ", i);
		matriceCovariance = covariances[i];
		for (j = 0; j < donnees->dimensionDesDonnees; j++) {
			for (k = 0; k < donnees->dimensionDesDonnees; k++) {
				fprintf(covariancesFile, "%f\t", matriceCovariance[j][k]);
			}
			fprintf(covariancesFile, "\n");
		}
		fprintf(covariancesFile, "\n");
	}
	fclose(covariancesFile);
	if (centres->dimensionDesDonnees > 2) {
		printf(
				"impossible de visualiser les données de cette dimension proceder à une reduction de dimension en ACP avant");
	} else {
		//generer un fichier .gnuplot et executer le avec system();
		FILE *gnuplotCmd = fopen("kmeans.gnu", "w");
		fprintf(gnuplotCmd,
				"set output 'kmeans.ps' \n set terminal postscript color\n plot 'centremobiles.txt' using 1:2 title 'centres', ");
		for (i = 0; i < config.centerNumber; i++) {
			char fileName[15];
			sprintf(fileName, "classe-%d.txt", i);
			fprintf(gnuplotCmd, "'%s' using 1:2 title 'classe-%d' ", fileName,
					i);
			if (i < (config.centerNumber - 1)) {
				fprintf(gnuplotCmd, ", ");
			}
			FILE *pointsParClasse = fopen(fileName, "w");
			for (j = 0; j < donnees->nombreDesDonnees; j++) {
				double* pointCourant = donnees->data[j];
				if (pointCourant[donnees->dimensionDesDonnees] == i) {
					for (k = 0; k < donnees->dimensionDesDonnees; k++) {
						fprintf(pointsParClasse, "%f\t", pointCourant[k]);
					}
					fprintf(pointsParClasse, "0\n");
				}
			}
			fprintf(pointsParClasse, "\n");
			fclose(pointsParClasse);
		}
		fprintf(gnuplotCmd, "\n");
		if (donnees->dimensionDesDonnees == 1) {
			fprintf(gnuplotCmd, "plot ");
			for (i = 0; i < config.centerNumber; i++) {
				fprintf(gnuplotCmd,
						"exp(-(x-%f)*(x-%f)/(2*%f))/sqrt(2*pi*%f) title 'gaussienne-%d' ",
						centres->codeBook[i][0], centres->codeBook[i][0],
						covariances[i][0][0], covariances[i][0][0], i);
				if (i < (config.centerNumber - 1)) {
					fprintf(gnuplotCmd, ", ");
				}
			}
		} else {
			//nombre de dimension = 2
			/*		fprintf(gnuplotCmd,"splot ");
			 for (i = 0; i < config.centerNumber; i++) {
			 matriceCovariance = covariances[i];
			 double sX = sqrt(matriceCovariance[0][0]);
			 double sY = sqrt(matriceCovariance[1][1]);
			 double covXY = matriceCovariance[0][1];
			 double coefficientDeCorrelation = covXY / (sX * sY);
			 fprintf(gnuplotCmd,
			 " exp(-1/2*(((x-%f)**2/%f)+((y-%f)**2/%f)-2*%f*((x-%f)*(y-%f)/(%f*%f))))/(2*pi*%f*%f*sqrt(1-%f**2)) title 'gaussienne-%d' ",
			 centres->codeBook[i][0], covariances[i][0][0],
			 centres->codeBook[i][1], covariances[i][1][1],
			 coefficientDeCorrelation, centres->codeBook[i][0],
			 centres->codeBook[i][1], sqrt(covariances[i][1][1]),
			 sqrt(covariances[i][0][0]), sqrt(covariances[i][1][1]),
			 sqrt(covariances[i][0][0]),coefficientDeCorrelation, i);
			 if (i < (config.centerNumber - 1)) {
			 fprintf(gnuplotCmd, ", ");
			 }
			 }
			 */}

		fprintf(gnuplotCmd,"\n set output\n");
		fclose(gnuplotCmd);

		system("gnuplot 'kmeans.gnu'");
	}
}
